NATRE region variance budget#
Reproducing the Ferrari & Polzin (2005) estimate using NATRE data, and comparing differences with our version
TODO#
[ ] Microscale variance budget
Could estimate \(⟨u_t T_t⟩ ⋅ ∇T_m = -⟨K_t T_z⟩ T_z^m\) term (monthly mean heat flux X monthly mean vertical gradient) using the model; just like I do with the microstructure data.
Then add the mesoscale dissipation term to get total \(⟨χ⟩\)
Setup#
Show code cell source
%load_ext watermark
%load_ext autoreload
%matplotlib inline
%matplotlib inline
import cf_xarray as cfxr # noqa
import cf_xarray.units # noqa
import colorcet
import dcpy
import distributed
import gsw
import holoviews as hv
import hvplot.xarray # noqa
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pint_xarray # noqa
import scipy as sp
import tqdm
hv.notebook_extension("bokeh")
import xgcm
from datatree import DataTree
from IPython.display import Image
import xarray as xr
from xarray.tests import raise_if_dask_computes
%aimport eddydiff
ed = eddydiff
xr.set_options(keep_attrs=True)
plt.rcParams["figure.dpi"] = 140
plt.rcParams["savefig.dpi"] = 200
plt.style.use("ggplot")
plt.figure()
plt.close()
%watermark -iv
xr.DataArray([1.0])
holoviews : 1.15.4
dcpy : 0.2.1.dev45+g67ec22a.d20230504
colorcet : 3.0.1
pint_xarray: 0.3
tqdm : 4.65.0
distributed: 2023.3.2
numpy : 1.23.5
hvplot : 0.8.3
matplotlib : 3.7.1
xarray : 2023.3.0
eddydiff : 0.1
gsw : 3.6.16
xgcm : 0.8.1
scipy : 1.10.1
cf_xarray : 0.8.1.dev0+gce9cd5a.d20230327
<xarray.DataArray (dim_0: 1)> array([1.]) Dimensions without coordinates: dim_0
Show code cell source
if "client" in locals():
client.cluster.close()
client.close()
client = distributed.Client(n_workers=3, threads_per_worker=1)
client
Client
Client-868c6dde-150a-11ee-ba4b-acde48001122
| Connection method: Cluster object | Cluster type: distributed.LocalCluster |
| Dashboard: http://127.0.0.1:8787/status |
Cluster Info
LocalCluster
e7a75469
| Dashboard: http://127.0.0.1:8787/status | Workers: 3 |
| Total threads: 3 | Total memory: 16.00 GiB |
| Status: running | Using processes: True |
Scheduler Info
Scheduler
Scheduler-5876b365-3b72-4551-9e25-7c9110c31f37
| Comm: tcp://127.0.0.1:54612 | Workers: 3 |
| Dashboard: http://127.0.0.1:8787/status | Total threads: 3 |
| Started: Just now | Total memory: 16.00 GiB |
Workers
Worker: 0
| Comm: tcp://127.0.0.1:54624 | Total threads: 1 |
| Dashboard: http://127.0.0.1:54628/status | Memory: 5.33 GiB |
| Nanny: tcp://127.0.0.1:54615 | |
| Local directory: /var/folders/yp/rzgd0f7n5zbcbf5dg73sfb8ntv3xxm/T/dask-worker-space/worker-oh8gs2vu | |
Worker: 1
| Comm: tcp://127.0.0.1:54623 | Total threads: 1 |
| Dashboard: http://127.0.0.1:54627/status | Memory: 5.33 GiB |
| Nanny: tcp://127.0.0.1:54616 | |
| Local directory: /var/folders/yp/rzgd0f7n5zbcbf5dg73sfb8ntv3xxm/T/dask-worker-space/worker-dx8t6h0s | |
Worker: 2
| Comm: tcp://127.0.0.1:54625 | Total threads: 1 |
| Dashboard: http://127.0.0.1:54626/status | Memory: 5.33 GiB |
| Nanny: tcp://127.0.0.1:54617 | |
| Local directory: /var/folders/yp/rzgd0f7n5zbcbf5dg73sfb8ntv3xxm/T/dask-worker-space/worker-lhjk4wik | |
Image("../images/natre-large-scale.png")
Prepare various production and dissipation terms#
Micro#
a05_binned = xr.load_dataset("../datasets/a05-2015.nc")
micro = xr.load_dataset("../datasets/natre-density-binned.nc")
Show code cell source
natre_fine = xr.load_dataset("../datasets/natre-finescale.nc")
natre_fine_avg = (
natre_fine.groupby_bins("pressure", bins=np.arange(250, 2001, 200))
.mean()
.cf.guess_coord_axis()
)
natre_fine_avg.pressure_bins.attrs = {
"positive": "down",
"standard_name": "sea_water_pressure",
}
Argo#
Show code cell source
argo_fine_profiles = xr.load_dataset(
"../datasets/argomix_natre.zarr", engine="zarr"
).drop("DIRECTION")
argo_fine_profiles["Kt"] = argo_fine_profiles["Kρ"]
argo_fine_profiles["Kt"].attrs["long_name"] = "$K_T$"
# argo_fine_profiles = argo_fine_profiles.where(argo_fine_profiles.ε < 0.1)
bins = cfxr.bounds_to_vertices(micro.gamma_n_bounds, "bound").data
argo_fine_dens = (
argo_fine_profiles.rename({"γmean": "gamma_n"})
.reset_coords("pressure")
.groupby_bins("gamma_n", bins=bins)
.mean()
)
argo_fine_dens["gamma_n_bins"].attrs = {
"standard_name": "neutral_density",
"axis": "Z",
"positive": "down",
}
argo_fine_dens["npts"] = (
argo_fine_profiles.χ.rename({"γmean": "gamma_n"})
.groupby_bins("gamma_n", bins=bins)
.count()
)
argo_fine_dens = argo_fine_dens.set_coords("pressure")
argo_fine_depth = argo_fine_profiles.groupby_bins(
"pressure", np.arange(250, 2050, 200)
).mean()
argo_fine_depth["pressure_bins"].attrs = {"standard_name": "sea_water_pressure"}
ECCO#
ecco_tile2 = xr.open_zarr("../datasets/ecco-tile2-sigma2.zarr")
ecco_tile2["THETA"].attrs["coordinates"] = "XC YC z_σ"
ecco_tile2["SALT"].attrs["coordinates"] = "XC YC z_σ"
ecco_natre = ecco_tile2.isel(i=slice(10, 18), j=slice(10, 20)).load()
ecco_natre_profile = ed.pop.calc_mean_redivar_profile(ecco_natre.mean("time"))
ecco_natre_initial = ed.pop.calc_mean_redivar_profile(ecco_natre.isel(time=0))
ecco_natre_final = ed.pop.calc_mean_redivar_profile(ecco_natre.isel(time=-1))
plt.legend(["mean", "initial", "final"])
ecco_natre_profile
<xarray.Dataset>
Dimensions: (σ: 60, bounds: 2)
Coordinates:
tile int64 2
z_σ (σ) float32 -64.37 -64.77 -65.33 ... -4.689e+03 -4.689e+03
* σ (σ) float32 34.15 34.15 34.17 34.18 ... 37.2 37.2 37.2 37.2
z_σ_bounds (σ, bounds) float32 -64.18 -64.57 ... -4.689e+03 -4.689e+03
Dimensions without coordinates: bounds
Data variables:
SALT (σ) float32 36.94 36.94 36.94 36.93 36.93 ... nan nan nan nan
THETA (σ) float32 21.72 21.69 21.65 21.6 21.47 ... nan nan nan nan
delT2_plane (σ) float64 1.35e-12 1.339e-12 1.326e-12 ... 0.0 0.0 0.0sub = ecco_tile2.SALT.sel(σ=36, method="nearest").load()
sub.sel(i=slice(40), j=slice(40)).hvplot.quadmesh(
geo=True,
widget_type="scrubber",
clim=(35.3, 35.8),
cmap="bmy_r",
widget_location="bottom",
)
ecco_natre_initial.delT2_plane.plot()
[<matplotlib.lines.Line2D at 0x17473b940>]
ecco = (
xr.open_dataset("../datasets/ecco-chi.nc")
.sel(time=slice("2000", None))
.mean("time")
)
for var in ecco.cf.coordinates["vertical"]:
ecco[var] = -1 * ecco[var]
kredi = ed.format_ecco2(
xr.open_dataset("~/datasets/ecco/interp_Ks_ECCOv4.nc")
).Ks_ECCOv4.rename({"pres": "k"})
ecco["Ke"] = (
kredi.sel(lat=slice(22, 28), lon=slice(322, 336)).mean(["lat", "lon"]).variable
)
ecco_natre_profile.update(ecco)
ecco_natre_profile["z_σ"] = -1 * ecco_natre_profile.z_σ
del ecco_natre_profile["σ"].attrs["positive"]
ecco_natre_profile["RediVar"] = ecco_natre_profile.GM_CHI
ecco_natre_profile
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/indexing.py:661: RuntimeWarning: deallocating CachingFileManager(<class 'netCDF4._netCDF4.Dataset'>, '/Users/dcherian/work/eddydiff/datasets/ecco-chi.nc', mode='r', kwargs={'clobber': True, 'diskless': False, 'persist': False, 'format': 'NETCDF4'}, manager_id='85ce065b-6ab5-4525-ac28-ceecad6c8c97'), but file is not already closed. This may indicate a bug.
self.array = NumpyIndexingAdapter(np.asarray(self.array))
<xarray.Dataset>
Dimensions: (σ: 60, node: 2, k: 50, bounds: 2, k_u: 50, k_l: 50, k_p1: 51)
Coordinates: (12/20)
tile int64 2
z_σ (σ) float32 64.37 64.77 65.33 ... 4.689e+03 4.689e+03 4.689e+03
* σ (σ) float32 34.15 34.15 34.17 34.18 ... 37.2 37.2 37.2 37.2
z_σ_bounds (σ, bounds) float32 -64.18 -64.57 ... -4.689e+03 -4.689e+03
* k (k) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
* k_u (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 41 42 43 44 45 46 47 48 49
... ...
drF (k) float32 10.0 10.0 10.0 10.0 ... 387.5 410.5 433.5 456.5
PHrefC (k) float32 49.05 147.1 245.2 ... 4.944e+04 5.357e+04 5.794e+04
PHrefF (k_p1) float32 0.0 98.1 196.2 ... 5.145e+04 5.57e+04 6.018e+04
rhoRef (k) float32 1.024e+03 1.024e+03 ... 1.052e+03 1.054e+03
drC_kl (k_l) float32 5.0 10.0 10.0 10.0 ... 376.0 399.0 422.0 445.0
* node (node) object 'adj' 'const'
Dimensions without coordinates: bounds
Data variables: (12/23)
SALT (σ) float32 36.94 36.94 36.94 36.93 36.93 ... nan nan nan nan
THETA (node, k) float32 22.61 22.59 22.54 22.38 ... 1.021 0.2023 0.0
delT2_plane (σ) float64 1.35e-12 1.339e-12 1.326e-12 ... 0.0 0.0 0.0
GM_CHIXO (node, k) float32 -1.491e-09 -2.558e-09 ... -7.767e-14 0.0
GM_CHIZ (node, k_l) float32 0.0 4.938e-09 5.708e-09 ... 1.499e-14 0.0
GM_ubT (node, k) float32 4.472e+05 -5.016e+04 ... 6.831e+03 0.0
... ...
Tx (node, k) float32 -1.636e-06 -1.639e-06 ... -1.105e-06 0.0
Ty (node, k) float32 -1.341e-06 -1.341e-06 ... -1.417e-08 0.0
Tz (node, k_l) float32 nan 0.001619 ... 0.001941 0.0004545
GM_CHI (node, k) float32 1.014e-08 9.531e-09 ... 4.292e-08 nan
Ke (k) float32 1.85e+03 2.002e+03 2.06e+03 ... 984.2 958.3 948.3
RediVar (node, k) float32 1.014e-08 9.531e-09 ... 4.292e-08 nanPOP#
POP 1/10° Variance budget#
Show code cell source
pop_hires_natre = (
xr.load_dataset("../datasets/natre-pop-variance.nc")
.pint.quantify()
.pint.to("K**2/s")
)
pop_hires_natre
<xarray.Dataset>
Dimensions: (depth: 62)
Coordinates:
* depth (depth) float64 5.0 15.0 25.0 ... 5.375e+03 5.625e+03 5.875e+03
Data variables:
BC (depth) float64 [K²/s] 3.765e-09 2.352e-09 ... -3.835e-12
HDIFF (depth) float64 [K²/s] -3.8e-10 -3.846e-10 ... -3.298e-14
PKC (depth) float64 [K²/s] 4.486e-10 6.327e-11 ... 4.232e-15 0.0
VMIX (depth) float64 [K²/s] -9.317e-09 -1.392e-09 ... 4.307e-14
DISS (depth) float64 [K²/s] -9.697e-09 -1.777e-09 ... 1.01e-14pop_hires_natre[["HDIFF", "VMIX"]].to_array().plot.line(
hue="variable", y="depth", ylim=(3000, 0)
)
[<matplotlib.lines.Line2D at 0x178a4ddf0>,
<matplotlib.lines.Line2D at 0x178a4de50>]
POP 1/10° state vars climatology#
pop_hires_clim = (
xr.load_dataset("../datasets/pop-hires-natre-coarsened-annual-climatology.nc")
.pint.quantify()
.squeeze("cycle")
.map(dcpy.util.to_base_units)
)
# pop_hires_clim['z_σ'] = pop_hires_clim.z_σ.cf.ffill("Z")
# pop_hires_clim = pop_hires_clim.cf.add_bounds("z_σ", dim="Z").isel(σ=slice(-7))
pop_hires_clim["TEMP"].attrs["standard_name"] = "sea_water_potential_temperature"
# NOTE: This only adds delT2, and averages in box
pop_hires_clim_profile = ed.pop.calc_mean_redivar_profile(pop_hires_clim)
# interpolate to pop_natre Z grid so we can estimate Ke
pop_hires_clim_profile = (
pop_hires_clim_profile.swap_dims({"σ": "z_σ"})
.drop_vars("z_σ_bounds")
.dropna("z_σ")
.interp(z_σ=pop_hires_natre.depth)
.cf.add_bounds("depth")
)
pop_hires_clim_profile["Ke"] = (
-1 * pop_hires_natre.DISS
) / pop_hires_clim_profile.delT2_plane.where(pop_hires_clim_profile.delT2_plane > 5e-13)
pop_hires_clim_profile
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/duck_array_ops.py:675: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return push(array, n, axis)
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/scipy/interpolate/_polyint.py:120: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
yi = np.asarray(yi)
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/scipy/interpolate/_interpolate.py:508: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
y = array(y, copy=self.copy)
<xarray.Dataset>
Dimensions: (depth: 62, bounds: 2)
Coordinates:
time object 0052-01-01 05:17:48.750000
σ (depth) float64 nan nan nan nan nan ... 37.2 nan nan nan nan
z_σ (depth) float64 5.0 15.0 25.0 ... 5.625e+03 5.875e+03
* depth (depth) float64 5.0 15.0 25.0 ... 5.625e+03 5.875e+03
depth_bounds (depth, bounds) float64 0.0 10.0 10.0 ... 5.75e+03 6e+03
z_σ_bounds (depth, bounds) float64 0.0 10.0 10.0 ... 5.75e+03 6e+03
Dimensions without coordinates: bounds
Data variables:
TEMP (depth) float64 nan nan nan nan nan ... 275.0 nan nan nan nan
SALT (depth) float64 nan nan nan nan nan ... nan nan nan nan
delT2 (depth) float64 nan nan nan nan nan ... nan nan nan nan
delT2_plane (depth) float64 nan nan nan nan nan ... nan nan nan nan
Ke (depth) float64 [K²/s] nan nan nan nan nan ... nan nan nan nan
Attributes: (12/14)
history: Mon Apr 15 11:06:38 2019: ncks -A 1st_half/ta...
title: g.e20.G.TL319_t13.control.001_hfreq
model_doi_url: https://doi.org/10.5065/D67H1H0V
time_period_freq: day_5
Conventions: CF-1.0; http://www.cgd.ucar.edu/cms/eaton/net...
yrs_averaged: 42-61
... ...
calendar: All years have exactly 365 days.
start_time: This dataset was created on 2019-01-16 at 20:...
contents: Diagnostic and Prognostic Variables
revision: $Id: tavg.F90 89091 2018-04-30 15:58:32Z altu...
history_of_appended_files: Mon Apr 15 11:06:38 2019: Appended file 1st_h...
NCO: netCDF Operators version 4.7.4 (http://nco.sf...POP 1°#
pop_1deg_natre_profile = xr.load_dataset(
"../datasets/pop-1deg-redivar-natre.nc"
).pint.quantify()
pop_1deg_natre_profile.σ.attrs["axis"] = "Z"
pop_1deg_natre_profile["z_σ"] = pop_1deg_natre_profile.z_σ.cf.ffill("Z")
pop_1deg_natre_profile = pop_1deg_natre_profile.cf.add_bounds("z_σ", dim="Z")
pop_1deg_natre_profile["delT2_plane"] = pop_1deg_natre_profile.delT2
pop_1deg_natre_profile
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/duck_array_ops.py:675: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return push(array, n, axis)
<xarray.Dataset>
Dimensions: (cycle: 6, σ: 60, bounds: 2)
Coordinates:
* cycle (cycle) int64 0 1 2 3 4 5
* σ (σ) float32 34.15 34.16 34.17 34.18 ... 37.2 37.2 37.2 37.2
z_σ (cycle, σ) float64 66.9 67.23 67.65 ... 5.036e+03 5.166e+03
z_σ_bounds (cycle, σ, bounds) float64 66.73 67.06 ... 5.101e+03 5.232e+03
Dimensions without coordinates: bounds
Data variables:
KAPPA_ISOP (cycle, σ) float64 [m²/s] 2.488e+03 2.486e+03 ... 508.3 464.9
RediVar (cycle, σ) float64 [K²/s] 2.956e-09 2.951e-09 ... 8.636e-13
TEMP (cycle, σ) float64 [K] 295.8 295.8 295.7 ... 274.8 274.8 274.8
delT2 (cycle, σ) float64 [K²/m²] 1.157e-12 1.156e-12 ... 1.605e-15
delT2_plane (cycle, σ) float64 [K²/m²] 1.157e-12 1.156e-12 ... 1.605e-15
Attributes:
Conventions: CF-1.0; http://www.cgd.ucar.edu/cms/eaton/netcdf/CF-cu...
calendar: All years have exactly 365 days.
cell_methods: cell_methods = time: mean ==> the variable values are ...
contents: Diagnostic and Prognostic Variables
history: none
model_doi_url: https://doi.org/10.5065/D67H1H0V
revision: $Id$
source: CCSM POP2, the CCSM Ocean Component
start_time: This dataset was created on 2019-06-20 at 13:53:17.7
time_period_freq: month_1
title: g.e21.GOMIPECOIAF_JRA.TL319_g17.CMIP6-omip2.001POP 1° spinup#
%autoreload
ds_ = xr.load_dataset(
"../datasets/pop-1deg-redi-var-natre-cycle-0.zarr/",
engine="zarr",
)
ds_.z_σ.attrs["units"] = "cm"
ds_.σ.attrs.update({"axis": "Z", "positive": "down"})
ds_.TEMP.attrs["standard_name"] = "sea_water_potential_temperature"
month1 = ds_.isel(yearmonth=0).pint.quantify().map(dcpy.util.to_base_units).squeeze()
month1_profile = ed.pop.calc_mean_redivar_profile(month1)
month1_profile = month1_profile.where(month1_profile.z_σ.notnull(), drop=True).isel(
σ=slice(1, None)
)
pop_1deg_spinup = (
ds_.squeeze("cycle")
.coarsen(yearmonth=120, boundary="trim")
.mean()
.rename({"yearmonth": "decade"})
.pint.quantify()
.map(dcpy.util.to_base_units)
)
pop_1deg_spinup_profile = ed.pop.calc_mean_redivar_profile(pop_1deg_spinup)
pop_1deg_spinup_profile
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/duck_array_ops.py:675: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return push(array, n, axis)
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/duck_array_ops.py:675: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return push(array, n, axis)
<xarray.Dataset>
Dimensions: (decade: 6, σ: 60, bounds: 2)
Coordinates:
cycle int64 0
month (decade) float64 6.5 6.5 6.5 6.5 6.5 6.5
time (decade) object 0005-12-16 00:00:00 ... 0055-12-16 00:00:00
year (decade) float64 5.5 15.5 25.5 35.5 45.5 55.5
* σ (σ) float32 34.15 34.15 34.17 34.18 ... 37.2 37.2 37.2 37.2
z_σ (decade, σ) float64 66.59 66.93 67.3 ... 4.931e+03 4.931e+03
z_σ_bounds (decade, σ, bounds) float64 66.42 66.76 ... 4.931e+03 4.931e+03
Dimensions without coordinates: decade, bounds
Data variables:
KAPPA_ISOP (decade, σ) float64 [m²/s] 2.512e+03 2.518e+03 ... nan nan
RediVar (decade, σ) float64 [K²/s] 3.03e-09 3.025e-09 ... nan nan
TEMP (decade, σ) float64 [K] 295.8 295.8 295.7 295.7 ... nan nan nan
delT2 (decade, σ) float64 [K²/m²] 1.186e-12 1.183e-12 ... nan nan
delT2_plane (decade, σ) float64 4.255e-13 4.351e-13 4.502e-13 ... 0.0 0.0
Attributes:
Conventions: CF-1.0; http://www.cgd.ucar.edu/cms/eaton/netcdf/CF-cu...
calendar: All years have exactly 365 days.
cell_methods: cell_methods = time: mean ==> the variable values are ...
contents: Diagnostic and Prognostic Variables
history: none
model_doi_url: https://doi.org/10.5065/D67H1H0V
revision: $Id$
source: CCSM POP2, the CCSM Ocean Component
start_time: This dataset was created on 2019-06-20 at 13:53:17.7
time_period_freq: month_1
title: g.e21.GOMIPECOIAF_JRA.TL319_g17.CMIP6-omip2.001Make DataTree#
pop_natre = DataTree()
pop_natre["1deg/cycle"] = DataTree(pop_1deg_natre_profile)
pop_natre["hires/clim"] = DataTree(
xr.merge([pop_hires_natre, pop_hires_clim_profile], join="exact")
)
pop_natre["1deg/decadal"] = DataTree(pop_1deg_spinup_profile)
pop_natre["1deg/month1"] = DataTree(month1_profile)
pop_natre
<xarray.DatasetView>
Dimensions: ()
Data variables:
*empty*Cole/Groeskamp#
This uses a crappier gradient.
Show code cell source
cole = ed.read_cole()
# read, take mean in region; swap to pressure coordinate
argograd = (
xr.open_dataset("../datasets/cole-clim-gradient-natre.nc")
.reset_coords("pres")
.sel(lat=slice(24, 28), lon=slice(360 - 30, 360 - 25))
.cf.guess_coord_axis()
)
argograd.lon.attrs["axis"] = "X"
argograd["lon"] = argograd.lon.cf.guess_coord_axis()
counts = argograd.CT.count(["lat", "lon"])
argo_delT2_plane = ed.plane_fit_gradient(argograd.CT).assign_coords(
{"pres": argograd.pres.mean(["lon", "lat"])}
)
cole_natre = cole[["diffusivity", "density_mean_depth", "depth_mean_sig"]].interp(
lat=argograd.lat.data, lon=argograd.lon.data
)
# move from sigma to depth space
argo_dTiso = xgcm.transform.linear_interpolation(
argograd.dTiso,
argograd.pres,
cole.pres,
"sigma",
"sigma",
"pres",
)
cole_natre["RediVar"] = cole_natre.diffusivity * argo_dTiso**2
cole_natre["RediVar"].attrs = {
"long_name": "$K_e^{cole} |∇T_{iso}^{argo}|²$",
"units": "°C²/s",
}
cole_natre = cole_natre.mean(["lat", "lon"])
groeskamp = (
xr.open_dataset("~/datasets/groeskamp2020/groeskamp2020.nc")
.sel(lat=slice(24, 28), lon=slice(360 - 30, 360 - 25))
.rename({"depth": "pres"})
)
groeskamp["dTiso"] = argo_dTiso.interp(pres=groeskamp.pres)
groeskamp["eddy_var_0"] = groeskamp["Ke_0"] * groeskamp.dTiso**2
groeskamp["eddy_var"] = groeskamp["Ke"] * groeskamp.dTiso**2
groeskamp = groeskamp.mean(["lat", "lon"])
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/indexing.py:661: RuntimeWarning: deallocating CachingFileManager(<class 'netCDF4._netCDF4.Dataset'>, '/Users/dcherian/datasets/groeskamp2020/groeskamp2020.nc', mode='r', kwargs={'clobber': True, 'diskless': False, 'persist': False, 'format': 'NETCDF4'}, manager_id='641c78cd-86a0-45ef-90b8-0a944d6638e8'), but file is not already closed. This may indicate a bug.
self.array = NumpyIndexingAdapter(np.asarray(self.array))
a05 = xr.load_dataset("../datasets/a05-2015.nc")
a05
<xarray.Dataset>
Dimensions: (bound: 2, gamma_n: 13, cast: 46)
Coordinates: (12/13)
* bound (bound) object 'lower' 'upper'
* gamma_n (gamma_n) float64 26.6 26.88 27.11 ... 27.93 27.96 27.98
num_obs (gamma_n) int64 3120 3281 3326 3256 ... 3549 3740 3463
pres (gamma_n) float64 233.1 383.6 ... 1.976e+03 2.149e+03
gamma_n_bounds (bound, gamma_n) float64 26.45 26.75 27.0 ... 27.97 27.99
station (cast) int64 108 109 110 111 112 ... 126 127 128 129 130
... ...
cleaner <U2 'GE'
time (cast) datetime64[ns] 2016-01-09T06:22:19.077246976 ......
lon (cast) float64 325.1 325.7 326.4 ... 337.7 338.4 339.0
lat (cast) float64 24.5 24.5 24.5 24.5 ... 24.92 25.13 25.33
CTD_chipod float64 9.969e+36
* cast (cast) int64 0 1 2 3 4 5 6 7 8 ... 38 39 40 41 42 43 44 45
Data variables: (12/59)
chi (gamma_n) float64 1.788e-08 7.233e-09 ... 2.8e-10
eps (gamma_n) float64 2.916e-09 1.259e-09 ... 2.938e-10
KtTz (gamma_n) float64 2.282e-07 2.495e-07 ... 6.227e-08
KtTz~ (gamma_n) float64 3.673e-07 1.921e-07 ... 6.11e-08
chib2 (gamma_n) float64 8.941e-09 3.617e-09 ... 1.4e-10
hm (gamma_n) float64 148.4 149.3 151.3 ... 165.5 173.1 162.3
... ...
δKtTzTz (gamma_n) float64 2.53e-09 2.664e-09 ... 5.358e-11
δKtTz~Tz (gamma_n) float64 1.783e-09 9.228e-10 ... 7.829e-11
δwTTz (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
δresidual (gamma_n) float64 -1.123e-09 -3.266e-10 ... -1.414e-11
δresidual_chi (gamma_n) float64 4.542e-10 1.32e-10 ... 2.716e-11
δpres (gamma_n) float64 74.2 74.67 75.67 ... 82.74 86.57 81.13
Attributes: (12/15)
cruise_name: A05
featureType: trajectoryProfile
Conventions: CF-1.8
title: Turbulence quantities measured by CTD-chipods mounted ...
summary: Turbulence data across the ocean basin from 10 m to th...
processing_level: 3: Variables mapped on uniform space-time grid from hi...
... ...
creator_name: Jonathan Nash / Aurelie Moulin
institution: Ocean Mixing Group, College of Earth, Ocean, and Atmos...
project: A05, GO-SHIP (The Global Ocean Ship-Based Hydrographic...
reference1: J.D. Nash et al., Ocean mixing measured by fast-respon...
comment: Expocode: 740H20200119
commit: eddydiff: b55f9526739dd63a768df2f816bc8327bcc0d247 | ...MIMOC#
mimoc = (
dcpy.oceans.read_mimoc(globstr="MIMOC_SIGMA0_*")
.coarsen(latitude=2, longitude=2, boundary="trim")
.mean()
)
mimoc.coords["lat_face"] = (
"lat_face",
mimoc.latitude.data - 0.25,
{"axis": "Y", "c_grid_axis_shift": -0.5},
)
mimoc.coords["lon_face"] = (
"lon_face",
mimoc.longitude.data - 0.25,
{"axis": "X", "c_grid_axis_shift": -0.5},
)
dx = gsw.distance(mimoc.lon_face.data, np.zeros_like(mimoc.lon_face.data))
dx = np.append(dx, dx[-1])
mimoc.coords["dx"] = ("longitude", dx, {"units": "m"})
mimoc.coords["dx"] = mimoc.dx * np.cos(np.pi / 180 * mimoc.latitude)
dy = gsw.distance(np.zeros_like(mimoc.latitude.data), mimoc.lat_face.data)
dy = np.append(dy, dy[-1])
mimoc.coords["dy"] = ("latitude", dy, {"units": "m"})
mimoc.coords["X"] = mimoc.dx.cf.cumsum("X")
mimoc.coords["Y"] = mimoc.dy.cf.cumsum("Y")
grid = xgcm.Grid(mimoc, metrics={"X": ["dx"], "Y": ["dy"]}, periodic=["X", "Y"])
dTdx = grid.derivative(mimoc.CT, "X")
dTdy = grid.derivative(mimoc.CT, "Y")
mimoc["delT2"] = grid.interp(dTdx, "X") ** 2 + grid.interp(dTdy, "Y") ** 2
mimoc.delT2.attrs = {"long_name": "$|∇T|^2$", "units": "C/m"}
mimoc_natre = mimoc.sel(
latitude=slice(25, 30), longitude=slice(360 - 30, 360 - 25)
).load()
mimoc_natre_profile = (
mimoc_natre.reset_coords("pressure")
.cf.mean(["X", "Y", "time"])
.set_coords("pressure")
)
mimoc_natre.sigma0.attrs["axis"] = "Z"
mimoc
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xgcm/grid.py:1645: UserWarning: Metric at ('time', 'lon_face', 'latitude', 'sigma0') being interpolated from metrics at dimensions ('longitude', 'latitude'). Boundary value set to 'extend'.
warnings.warn(
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xgcm/grid.py:1645: UserWarning: Metric at ('time', 'longitude', 'lat_face', 'sigma0') being interpolated from metrics at dimensions ('latitude',). Boundary value set to 'extend'.
warnings.warn(
<xarray.Dataset>
Dimensions: (time: 12, longitude: 360, latitude: 170, sigma0: 75,
lat_face: 170, lon_face: 360)
Coordinates:
* latitude (latitude) float32 -79.75 -78.75 -77.75 ... 87.25 88.25 89.25
* longitude (longitude) float32 0.25 1.25 2.25 3.25 ... 357.2 358.2 359.2
* sigma0 (sigma0) float32 22.4 22.64 22.88 23.11 ... 28.14 28.15 28.24
* time (time) datetime64[ns] 2014-01-15 2014-02-15 ... 2014-12-15
* lat_face (lat_face) float32 -80.0 -79.0 -78.0 -77.0 ... 87.0 88.0 89.0
* lon_face (lon_face) float32 0.0 1.0 2.0 3.0 ... 357.0 358.0 359.0
dx (longitude, latitude) float64 1.979e+04 ... 1.455e+03
dy (latitude) float64 1.112e+05 1.112e+05 ... 1.112e+05
X (longitude, latitude) float64 1.979e+04 ... 5.24e+05
Y (latitude) float64 1.112e+05 2.224e+05 ... 1.879e+07 1.89e+07
Data variables:
SUMMED_WEIGHT (time, longitude, latitude, sigma0) float32 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
YEAR_OF_DATA (time, longitude, latitude, sigma0) float32 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
pressure (time, longitude, latitude, sigma0) float32 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
SA (time, longitude, latitude, sigma0) float32 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
CT (time, longitude, latitude, sigma0) float32 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
delT2 (time, longitude, latitude, sigma0) float64 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>mimoc_natre.CT.isel(sigma0=24, time=0).cf.plot()
<matplotlib.collections.QuadMesh at 0x178bd8e20>
mimoc_natre_profile["delT2_plane"] = ed.plane_fit_gradient(
mimoc_natre.CT.mean("time"), debug=True
)
argo_delT2_plane = ed.plane_fit_gradient(argoclim.CT, debug=True).assign_coords(
{"pres": argoclim.pres.mean(["lon", "lat"])}
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 4
1 mimoc_natre_profile["delT2_plane"] = ed.plane_fit_gradient(
2 mimoc_natre.CT.mean("time"), debug=True
3 )
----> 4 argo_delT2_plane = ed.plane_fit_gradient(argoclim.CT, debug=True).assign_coords(
5 {"pres": argoclim.pres.mean(["lon", "lat"])}
6 )
NameError: name 'argoclim' is not defined
AVISO#
aviso = xr.open_zarr("~/datasets/aviso_monthly.zarr", consolidated=True)
aviso["longitude"] = aviso.longitude - 360
aviso = aviso.sel(latitude=slice(10, 50), longitude=slice(-40, 5))
aviso
Variance production vs dissipation#
Agreement between POP “mesoscale variance dissipation” and microstructure estimate of “variance received from mesoscale” is great! (Purple dashed steps agree well with purple solid line)
POP’s grid-scale temp variance dissipation is really dominated by the horizontal biharmonic diffusion, but fundamentally the amount of mesoscale variance destroyed by grid-scale diffusion matches that in the microstructure measurements.
One place of disagreement is between 400m and 700m where the POP estimate is really small compared to my difference. FP2005 also see this using the mooring based Ke estimate. At these levels their two terms also match a lot closer than mine, so there’s an error somewhere. It seems like my \(⟨χ⟩/2\) is a little bit too high which happened after I switched to the bootstrapping based estimates. I’ve made some improvements
Sorting the dataset, changes the mean density field used to pick the bins.
Calculating \((∂T/∂z)\_m\) is a little wonky. I now use a different procedure. In each O(100m) \(γ_n\) bin, I groupby finer-resolution \(γ_n\) bins (21 of them), calculate the mean profile, then fit a straight line.
Cole/Groeskamp \(K_e\) estimate based variance dissipation compares much better now.
~Comparing to variance dissipated by parameterizations for \(K_e\) (\(K_e |∇_hT|²\)) doesn’t work well. I wonder if I’m doing something wrong with my \(∇_hT\) estimate~
I updated my |∇T| estimation code; now estimate along-isopycnal gradients in σ space instead of just a lateral gradient
While the current algo still doesn’t match Silvia’s ∇S (but much better than earlier), it seems to improve agreement (panel c)
Still a lot of disagreement below 1200m when using published \(K_e\) estimates.
Show code cell source
plt.rcParams["figure.dpi"] = 140
plt.rcParams["figure.facecolor"] = "none" # (0.91, 0.91, 0.91)
plt.rcParams["axes.facecolor"] = (0.99,) * 3
plt.rcParams.update(
{
"axes.spines.top": True,
"axes.spines.bottom": True,
"axes.spines.left": True,
"axes.spines.right": True,
"axes.edgecolor": (0.31,) * 3,
"axes.grid.which": "both",
"grid.linewidth": 0.5,
"grid.alpha": 0.3,
"grid.color": (0,) * 3,
}
)
micro.pres.attrs["bounds"] = "pres_err"
micro["residual"] = np.abs(micro.residual)
micro["residual_err"] = np.abs(micro.residual_err)
hires = pop_natre["hires/clim"].drop_vars(["z_σ", "z_σ_bounds"])
mesomicrocolor = "darkorchid"
error_kwargs = dict(error="x", ls="none", marker=".")
def plot_pop_natre_meso_micro(ax):
(-1 * hires["DISS"]).cf.plot(
ax=ax, lw=1.5, color=mesomicrocolor, label="POP meso $→$ micro"
)
dcpy.plots.errorbar(
micro,
x="residual",
y="pres",
color=mesomicrocolor,
label="NATRE: $⟨χ⟩/2 - K_ρ^m ∂_zθ_m^2$",
# lw=2.5,
ax=ax,
**error_kwargs,
)
# dcpy.plots.fill_between_bounds(
# micro,
# "residual",
# y="pres",
# color=mesomicrocolor,
# label="NATRE meso $→$ micro",
# ax=ax,
# )
# f, axx = plt.subplots(1, 4, )
f, ax = plt.subplot_mosaic(
[["natre", "pop-hires", "params"]],
sharey=False,
sharex=True,
constrained_layout=False,
)
### Mesoscale budget
hires["BC"].cf.plot(ax=ax["pop-hires"], label="POP h. stir $⟨u_e^h θ_e⟩.∇_hθ_m$")
np.abs(-1 * hires["PKC"]).cf.plot(
ax=ax["pop-hires"], label="POP v. stir $-⟨w_eθ_e⟩.∂_zθ_m$"
)
# (-1 * pop_natre.VMIX).cf.plot(color="sienna", ax=ax[0], label="POP VMIX")
# (-1 * pop_natre.HDIFF).cf.plot(ax=ax[0], label="POP HDIFF", marker="x", ls="none")
#################################
######## Microscale budget: NATRE
#################################
dcpy.plots.fill_between_bounds(
micro,
"KρTz2",
y="pres",
color="k",
label="FP2005: $K_ρ^m ∂_zθ_m^2$",
ax=ax["natre"],
)
# dcpy.plots.fill_between_bounds(
# micro,
# "KtTz~Tz",
# y="pres",
# color="C1",
# label="NATRE $⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
# ax=ax[1],
# )
dcpy.plots.fill_between_bounds(
micro,
"chib2",
y="pres",
color="C0",
label="FP2005: $⟨χ⟩/2$",
ax=ax["natre"],
)
# (argo_fine_depth.χ.sel(criteria="whalen") / 2).cf.plot.step(
# # marker="o",
# color="green",
# ax=ax[1],
# y="pressure_bins",
# label="$⟨χ⟩^{argo}_{Rω=3}/2$",
# )
#################################
########## Microscale budget: A05
#################################
if "a05" in ax:
dcpy.plots.fill_between_bounds(
a05_binned,
"KtTz~Tz",
y="pres",
color="k",
label=r"$⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
ax=ax["a05"],
)
dcpy.plots.fill_between_bounds(
a05_binned,
"chib2",
y="pres",
color="C0",
label="$⟨χ⟩/2$",
ax=ax["a05"],
)
dcpy.plots.errorbar(
a05_binned,
x="residual_chi",
y="pres",
color=mesomicrocolor,
label="A05 meso → micro",
ax=ax["a05"],
**error_kwargs,
)
#################################
### Parameterizations
#################################
cole_natre.RediVar.plot(
y="pres", label="$K_e^{cole} |∇_hT^{cole}|²$", color="k", ls="--", ax=ax["params"]
)
# dcpy.plots.errorbar(
# micro,
# x="residual",
# y="pres",
# # color=mesomicrocolor,
# label="NATRE $⟨χ⟩/2 - ⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
# ax=ax[3],
# **error_kwargs,
# )
if "a05" in ax:
dcpy.plots.errorbar(
a05_binned,
x="residual_chi",
y="pres",
# color=mesomicrocolor,
label=r"A05 $⟨χ⟩/2 - ⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
ax=ax["params"],
**error_kwargs,
)
# groeskamp.eddy_var_0.plot(
# y="pres",
# label="$K_e^{G2020unsupp} |∇_hT^{cole}|²$",
# color="limegreen",
# ls="-.",
# ax=ax[2],
# )
groeskamp.eddy_var.plot(
y="pres",
label="$K_e^{G2020} |∇_hT^{cole}|²$",
color="limegreen",
ls="--",
ax=ax["params"],
)
#################################
### Cleanup
#################################
for axx in ax.values():
plot_pop_natre_meso_micro(axx)
axx.set_xscale("log")
axx.set_xlabel("")
axx.set_ylim([1850, 350])
axx.set_xlim([1e-11, 1e-8])
# axx.set_xticks([1e-11, 1e-10, 1e-9, 1e-8])
axx.xaxis.set_minor_locator(
mpl.ticker.LogLocator(base=10, subs="all", numticks=120)
)
# axx.set_ylim([1900, 180])
axx.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15))
axx.grid(True, which="both", lw=0.5)
patch = mpl.patches.Rectangle(
(0, 800), 1, 700, transform=axx.get_yaxis_transform(), alpha=0.1, zorder=1
)
axx.add_patch(patch)
# ax[1].set_xlabel("Variance production or dissipation [°C²/s]")
# f.suptitle("NATRE region")
ax["pop-hires"].set_title("Mesoscale budget: POP2 1/10°", fontsize="medium")
ax["natre"].set_title("Microscale budget: NATRE", fontsize="medium")
if "a05" in ax:
ax["a05"].set_title("Microscale budget: A05", fontsize="medium")
ax["params"].set_title("$K_e$ Parameterizations", fontsize="medium")
dcpy.plots.label_subplots(ax.values(), x=0.03, y=0.95)
dcpy.plots.clean_axes(list(ax.values()))
f.set_size_inches((8, 5))
f.savefig("../images/natre-meso-micro-param.pdf", bbox_inches="tight")
f.savefig("../images/natre-meso-micro-param.png", bbox_inches="tight")
Show code cell source
micro.pres.attrs["bounds"] = "pres_err"
mesomicrocolor = "darkorchid"
error_kwargs = dict(error="x", ls="none", marker=".")
def plot_pop_natre_meso_micro(ax):
(-1 * pop_natre.DISS).cf.plot(
ax=ax, lw=1.5, color=mesomicrocolor, label="POP meso $→$ micro"
)
# dcpy.plots.fill_between_bounds(
# micro,
# "residual",
# y="pres",
# color=mesomicrocolor,
# label="NATRE meso $→$ micro",
# ax=ax,
# )
f, ax = plt.subplots(1, 3, sharey=False, sharex=True, constrained_layout=False)
### Mesoscale budget
# pop_natre.BC.cf.plot(ax=ax[0], label="POP h. stir $⟨u_e^h θ_e⟩.∇_hθ_m$")
# np.abs(-1 * pop_natre.PKC).cf.plot(ax=ax[0], label="POP v. stir $-⟨w_eθ_e⟩.∂_zθ_m$")
# (-1 * pop_natre.VMIX).cf.plot(color="sienna", ax=ax[0], label="POP VMIX")
# (-1 * pop_natre.HDIFF).cf.plot(ax=ax[0], label="POP HDIFF", marker="x", ls="none")
#################################
######## Microscale budget: NATRE
#################################
plt.sca(ax[0])
dcpy.plots.fill_between_bounds(
micro,
"KρTz2",
y="pres",
color="k",
label="FP2005: $χ_m = K_ρ^m ∂_zθ_m^2$",
ax=ax[0],
)
# dcpy.plots.fill_between_bounds(
# micro,
# "KtTz~Tz",
# y="pres",
# color="C1",
# label="NATRE $⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
# ax=ax[0],
# )
dcpy.plots.fill_between_bounds(
micro,
"chib2",
y="pres",
color="C0",
label="FP2005: $⟨χ⟩/2$",
ax=ax[0],
)
dcpy.plots.errorbar(
micro,
x="residual",
y="pres",
color=mesomicrocolor,
label="NATRE $χ_e$",
ax=ax[0],
**error_kwargs,
)
# (argo_fine_depth.χ.sel(criteria="whalen") / 2).cf.plot.step(
# # marker="o",
# color="green",
# ax=ax[1],
# y="pressure_bins",
# label="$⟨χ⟩^{argo}_{Rω=3}/2$",
# )
#################################
########## Microscale budget: A05
#################################
dcpy.plots.fill_between_bounds(
a05_binned,
"KtTz~Tz",
y="pres",
color="k",
label="$χ_m$",
ax=ax[1],
)
dcpy.plots.fill_between_bounds(
a05_binned,
"chib2",
y="pres",
color="C0",
label="$⟨χ⟩/2$",
ax=ax[1],
)
dcpy.plots.errorbar(
a05_binned,
x="residual_chi",
y="pres",
color=mesomicrocolor,
label="A05 $χ_e$",
ax=ax[1],
**error_kwargs,
)
#################################
### Parameterizations
#################################
cole_natre.RediVar.plot(
y="pres", label="$K_e^{cole} |∇_hT^{cole}|²$", color="k", ls="--", ax=ax[2]
)
dcpy.plots.errorbar(
micro,
x="residual",
y="pres",
# color=mesomicrocolor,
label="NATRE: $χ_e$",
lw=2.5,
ax=ax[2],
**error_kwargs,
)
# dcpy.plots.errorbar(
# micro,
# x="residual",
# y="pres",
# # color=mesomicrocolor,
# label="NATRE $⟨χ⟩/2 - ⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
# ax=ax[3],
# **error_kwargs,
# )
dcpy.plots.errorbar(
a05_binned,
x="residual_chi",
y="pres",
# color=mesomicrocolor,
label="A05 $χ_e$",
ax=ax[2],
**error_kwargs,
)
# groeskamp.eddy_var_0.plot(
# y="pres",
# label="$K_e^{G2020unsupp} |∇_hT^{cole}|²$",
# color="limegreen",
# ls="-.",
# ax=ax[2],
# )
groeskamp.eddy_var.plot(
y="pres", label="$K_e^{G2020} |∇_hT^{cole}|²$", color="limegreen", ls="--", ax=ax[2]
)
#################################
### Cleanup
#################################
for axx in ax:
plot_pop_natre_meso_micro(axx)
axx.set_xscale("log")
axx.set_xlabel("")
axx.set_ylim([1850, 350])
axx.set_xlim([1e-11, 1e-8])
# axx.set_xticks([1e-11, 1e-10, 1e-9, 1e-8])
axx.xaxis.set_minor_locator(
mpl.ticker.LogLocator(base=10, subs="all", numticks=120)
)
# axx.set_ylim([1900, 180])
axx.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15))
axx.grid(True, which="both", lw=0.5)
# ax[1].set_xlabel("Variance production or dissipation [°C²/s]")
f.suptitle("NATRE region")
# ax[0].set_title("(a) Mesoscale budget", fontsize="medium")
ax[0].set_title("(a) Microscale budget: NATRE", fontsize="medium")
ax[1].set_title("(b) Microscale budget: A05", fontsize="medium")
ax[2].set_title("(c) Parameterizations", fontsize="medium")
dcpy.plots.clean_axes(ax)
f.set_size_inches((8, 5))
f.savefig("../images/natre-3-panel.png", bbox_inches="tight")
Interpretation#
Our analytical framework follows that of Ferrari & Polzin (2005) and Garrett (1992). Begin with a triple decomposition of the temperature fields
where the subscripts \(m\), \(e\) and \(t\) are loosely identified with the “mean”, “mesoscale eddy”, and “turbulence” scales.
Define the “large-scale averaging operator” \(⟨⟩\) that filters out the “eddy” and “turbulence” so that \(⟨T⟩ = T_m\) and an intermediate scale averaging operator \(\widetilde{\phantom{...}}\) that filters out the “turbulence” so that \(\widetilde{T} = T_m + T_e\).
Ferrari & Polzin (2005) follow Garrett (1992) to derive the variance budget equations for \(⟨T_e^2⟩\) and \(⟨T_t^2⟩\). Assuming stationarity, and homogeneity, and /a priori/ ignoring the triple correlation term, we can write the approximate equations
Here \(χ\) representes ths instantaneous rate of dissipation of termpeature variance and \(κ_T\) is the molecular diffusiovity of temperature.
The first equation, for mesoscale variance \(⟨T_e^2⟩\), states that the mesoscale scale field generates variance by the stirring of the mean and this variance cascades down to the microscale turbulence through the “scale transoframtion term” \(⟨\widetilde{u_t T_t} ⋅ ∇T_e⟩\). The second equation, for microscale variance \(⟨T_t^2⟩\), states that the microscale turbulence stirs the combination of the mean and eddy fields \((T_m + T_e)\) to generate variance that is eventually dissipated at the moelcular scale at the rate \(χ\). The scale transformation term \(⟨\widetilde{u_t T_t} ⋅ ∇T_e⟩ \)has opposite sign in the two equations indicating that it is a link between the mesoscale and microscale. Adding the two equations yields the approximate blaance,
that is the variance dissipated at the molecular scale is approximately generated by the stirring of the mean field \(T_m\) by the mesoscale and microscale turbulence.
A large number of approximations are required to get to this point, but this framework has value in qualitatively describing the nature of variance cascades in the ocean \citep{Garrett1992,Davis1994,Ferrari2005,NaveriaGarabato2016,Springys2020}.
Spicy region#
We can interpret the red mesoscale dissipation term in (b) as the amount of variance transferred to the microscale for eventual dissipation. This peaks in regions of high mean spice between 1000m and 1500m. The balance in spicy regions is that horizontal mesoscale stirring converts MPE → EPE (compensated filaments or T-S variance) that is dissipated by the microscale. The approximate budget is (note vertical stirring gets ignored, and \(u^h\) represents horizontal velocity vector)
Yiming sees the first equation as “lateral stirring” balances “mesoscale dissipation”, and
FP2005 combine the two equations to conclude that “mesoscale stirring of mean variance” balances “total microscale dissipation”
Weak spice#
In the upper water column that is not spicy, we instead see horizontal stirring MPE→ EPE which is then routed to EKE by “vertical stirring”. Here the scale transformation term is 0 so the budget is
Yiming sees the first equation as “low dissipation” or “lateral stirring balances vertical stirring”;
FP2005 sees the second equation as “microscale stirring of the mean” balances “total microscale dissipation”
All microscale budget terms#
Microscale budget:
Ferrari & Polzin use \(⟨u_t θ_t⟩.∇θ_m = -Γ ⟨ε⟩/N_m^2 ∂_zθ_m^2\), we are instead doing \(-⟨K_T θ_z⟩ ∂_zθ_m\).
Our estimate mostly agrees with theirs (crosses versus solid-line-steps).
Filling in NaNs in the temperature profiles was important. This means more estimates are used in the averages
Masked where both \(χ\) and \(ε\) are both available.
Masked so that \(|θ_z|\) > 1e-3 (really important, otherwise our estimate is too high)
Variance produced by parameterizations for eddy diffusivity:
The comparison with variance production using Cole et al diffusivity and gradients calculated in a similar way disagree significantly between 400 and 1000m. Below that there’s not much data but it seems to agree.
The Groeskamp et al (2020) estimate [G2020] is slightly better but not by much. All predict a maximum in variance production near 800m
⟨χ⟩ estimate using finescale parameterizations:
The finescale estimates are interesting, argo with Rω=3 and natre with Rω=7 work well but why?
There is some interesting discrepancy here. This choice of \(Rω\) makes for a factor of 2 difference which would hide the mesoscale→microscale transformation.
f, ax = plt.subplots(1, 1) # , constrained_layout=True)
(micro.chi / 2).cf.plot.step(
y="pres",
color="r",
label="FP: $⟨χ⟩/2$",
lw=2,
xlim=(1e-11, 1e-8),
ylim=(1800, 200),
)
(micro.Krho_m * micro.dTdz_m**2).cf.plot.step(
color="k", lw=2, y="pres", label="FP: $Γ ⟨ε⟩/N_m^2 ∂_zθ_m^2$"
)
cole_natre.RediVar.plot(y="pres", label="$K_e^{cole} |∇T^{cole}|²$", color="k", ls="--")
groeskamp.eddy_var_0.plot(
y="depth", label="$K_e^{G2020unsupp} |∇T^{cole}|²$", color="teal", ls="-."
)
groeskamp.eddy_var.plot(
y="pres", label="$K_e^{G2020} |∇T^{cole}|²$", color="teal", ls="--"
)
ed.sections.plot_var_prod_diss(micro, ls="none", marker="x", ax=plt.gca())
(argo_fine_depth.χ.sel(criteria="whalen") / 2).cf.plot.step(
marker="o",
color="C1",
ax=plt.gca(),
y="pressure_bins",
label="$⟨χ⟩^{argo}_{Rω=3}/2$",
)
# (natre_fine_avg.χ.sel(criteria="whalen") / 2).cf.plot.step(
# marker="o", color="C2", ax=plt.gca(), label="$⟨χ⟩_{fine-3}^{natre}/2$"
# )
(natre_fine_avg.χ.sel(criteria="whalen_7") / 2).cf.plot.step(
marker="o", color="C3", ax=plt.gca(), label="$⟨χ⟩_{Rω=7}^{natre}/2$"
)
# (micro.chi / 2 - micro.Krho * micro.dTdz ** 2).cf.plot.step(
# y="pres",
# xscale="log",
# xlim=(1e-11, 1e-8),
# ylim=(1800, 200),
# color="darkblue",
# lw=1,
# label="$⟨χ⟩/2 - K_ρ θ^m_z$",
# )
plt.legend(loc="right", bbox_to_anchor=(1.6, 0.5))
plt.xlabel("Variance production or dissipation [°C²/s]")
plt.gcf().set_size_inches((4.5, 6))
plt.title("NATRE")
plt.gcf().savefig("../images/natre-estimate.png")
POP2 1/10° mesoscale variance budget#
Yiming is working with the following variance budget equation for the mesoscale after decomposing \(T = \bar{T} + T' + T'' ≡ T_m + T_e + T_t\) (the latter is Ferrari & Polzin notation). The \(T = \tilde{T} + T''\) decomposition is inherent when using cell-averaged equations.
Ferrari & Polzin’s equations (4.6, 4.7) for the mesoscale and microscale are (\(u\) here is the full 3D \(u\)):
Note that these equations are specialized for the “interior”.
For the model the averaging scales are:
~: grid-scale in vertical; single time-step in time i.e. the “microscale” is the “grid-scale”
⟨⟩: monthly in time; NATRE domain size in (x,y)
In Yiming’s notation these equations are
Now we can identify Yiming’s dissipation term (“VMIX” + “HDIFF”, ignoring solar heating etc.) as the scale transformation term in FP2005 and then use the microscale budget from FP2005 to write
i.e. Yiming’s “mesoscale dissipation” term is really the difference between
total microscale dissipation, and
microscale production by stirring of the mean vertical gradient,
both of which can be estimated from microstructure observations.
Show code cell source
Image("../images/Tvariance_budget_NATREregion_with_EKE.png", width=600)
distributed.utils_perf - WARNING - full garbage collections took 13% CPU time recently (threshold: 10%)
Here are Yiming’s equations
Image("../yiming-budget.png")
POP 1/10 vs POP 1 vs ECCO#
TODO:#
[ ] Add ECCO IC
[ ] interpolate diffuivities to density surfaces?
[ ] Check
pop-hires-annual-climatology.nc[ ] For 1/10°, Infer 1° K_e, then average…?
[ ] Add Scott’s K; diagonal components; maybe compare cross-stream diffusivities?
[ ] Some tuning of Kappa for SMYLE.
[ ] 1/10° BGC run from MATT
[x] Email Keith about initialization for 1°.
[x] Compare ∇T for POP 1° and 1/10° initial conditions.
This plot compare POP simulations.
POP 1/10° is regridded to 1° grid.
All model output is averaged over years 42-61
1° grid products are regridded to density space; then averaged in NATRE box for mean \(K_e, χ^{meso}\)
For \(|∇_ρT|\) I take 1° T field; and do plane fit in NATRE box; then estimate \(|∇T|^2\)
Conclusions:#
POP 1/10° |∇T| agrees decently well with Argo |∇T| ; Argo and MIMOC agree quite well
POP 1° initial conditions don’t capture the mean spice gradient in the NATRE region.
1/10° initialized with quarter degree WOA?
1° initialized with something else?
POP 1° \(K\_{redi}\) is too high, below a 1000m, by a factor of 2-3.
Inferring \(K_e = χ/|∇T|^2\) for the 1° gives basically the same value as
KAPPA_ISOPPurely diffusive view of the problem. What if the production - med outflow - is too weak?
argo_delT2_plane.pres.attrs["positive"] = "down"
Show code cell source
plt.rcParams["figure.dpi"] = 140
plt.rcParams["figure.facecolor"] = "none" # (0.91, 0.91, 0.91)
plt.rcParams["axes.facecolor"] = (0.99,) * 3
plt.rcParams.update(
{
"axes.spines.top": True,
"axes.spines.bottom": True,
"axes.spines.left": True,
"axes.spines.right": True,
"axes.edgecolor": (0.31,) * 3,
"axes.grid.which": "both",
"grid.linewidth": 0.5,
"grid.alpha": 0.3,
"grid.color": (0,) * 3,
}
)
hirescolor = "darkorchid"
spinupcolor = "r"
spinupendcolor = "b"
lastcyclecolor = "deepskyblue"
color_ecco = "darkorange"
f, axx = plt.subplots(1, 3, sharey=True, constrained_layout=True)
ax = dict(zip(["K", "delT2", "chi"], axx))
legend_handles = []
for pop_, color, label in zip(
[
pop_natre["hires/clim"],
# pop_1deg_spinup_profile.isel(decade=0),
pop_natre["1deg/month1"],
pop_natre["1deg/decadal"].isel(decade=5),
# pop_1deg_natre_profile.sel(cycle=5),
],
[hirescolor, spinupcolor, spinupendcolor, lastcyclecolor],
["POP2 1/10°", "POP2 1° month 1", "POP2 1° decade 6", "POP2 1° last cycle"],
):
kwargs = dict(
edges=ed.pop.get_edges(pop_.ds), orientation="horizontal", lw=2, color=color
)
ax["delT2"].stairs(np.sqrt(pop_["delT2_plane"].data), **kwargs)
if "KAPPA_ISOP" in pop_:
legend_handles.append(ax["K"].stairs(pop_["KAPPA_ISOP"].to_numpy(), **kwargs))
ax["chi"].stairs(pop_["RediVar"].to_numpy(), **kwargs, label=label)
# K
legend_handles.append(
groeskamp.Ke.cf.plot(
y="pres", ax=ax["K"], color="limegreen", ls="--", label="Groeskamp (2020)"
)
)
legend_handles.append(
cole_natre.diffusivity.cf.plot(
y="pres", ax=ax["K"], color="k", ls="--", label="Cole (2015)"
)
)
legend_handles.append(ecco.Ke.cf.plot(y="vertical", ax=ax["K"], color=color_ecco))
# pop_hires_clim_profile.Ke.drop_vars("z_σ").cf.plot.step(
# y="depth",
# ax=ax["K"],
# xlim=[0, 3000],
# color=hirescolor,
# label="POP 1/10° χ$^{e}$/|∇T|$^2$",
# )
# ∇T2
for data, color, label in zip(
[
argo_delT2_plane, # mimoc_natre_profile.delT2_plane.rename({"pressure": "pres"})
ecco_natre_profile.delT2_plane,
],
[
"k",
color_ecco, # "darkcyan"
],
[
"Argo",
"ECCO Adjusted", # "MIMOC"
],
):
if label == "MIMOC":
continue
np.sqrt(data).cf.plot(
y="vertical", ax=ax["delT2"], color=color, label=label, _labels=False
)
# χmeso
(-1 * pop_natre["hires/clim"]["DISS"].drop_vars(["σ", "z_σ"])).cf.plot(
ax=ax["chi"], lw=2, color=hirescolor, label="POP2 1/10°", _labels=False
)
(ecco.GM_CHI).sel(node="adj").cf.plot(
ax=ax["chi"],
y="vertical",
lw=2,
color=color_ecco,
label=None,
_labels=False,
)
# ecco.GM_CHI.sel(node="const").cf.plot(
# ax=ax["chi"], y="vertical", lw=2, label="ECCO const", _labels=False
# )
dcpy.plots.errorbar(
micro,
x="residual",
y="pres",
color=mesomicrocolor,
label="NATRE: $⟨χ⟩/2 - K_ρ^m ∂_zθ_m^2$",
# lw=2.5,
ax=ax["chi"],
**error_kwargs,
)
ax["K"].set_title("")
ax["K"].set_ylim([3000, 200])
ax["K"].set_xlim([0, 2000])
ax["delT2"].set_xlim([1e-8, 4e-6])
ax["chi"].set_xlim([1e-12, 1e-8])
dcpy.plots.clean_axes(axx)
for a in axx:
# a.legend(loc="upper center", bbox_to_anchor=(0.5, -0.0), frameon=False)
a.xaxis.tick_top()
a.xaxis.set_label_position("top")
patch = mpl.patches.Rectangle(
(0, 800), 1, 700, transform=a.get_yaxis_transform(), alpha=0.1, zorder=1
)
a.add_patch(patch)
ax["K"].set_xlabel("(a) $K_{e}$ [m$^2$/s]", labelpad=12)
ax["delT2"].set_xlabel("(b) $|∇_ρ^hT|^2$ [C/m]", labelpad=12)
ax["chi"].set_xlabel("(c) $χ_{e} = K_{e} |∇_ρ^hT|^2$ [C$^2$/s]", labelpad=12)
ax["K"].set_ylabel("depth [m]")
ax["chi"].set_xscale("log")
ax["chi"].xaxis.grid(True, which="both")
f.legend(loc="outside right")
# dcpy.plots.label_subplots(list(ax.values()), x=0.9, y=0.1)
f.set_size_inches((9, 6))
f.savefig("../images/coarse-models.pdf")
f.savefig("../images/coarse-models.png")
ECCO#
Kux, Kvy are constant in time as expected.
But when I turn off the adjusted file, I still get values ≠ 1000
Interestingly adjusted \(K\) is higher right where we expect variance dissipation.
Adjust \(K\) is also high right where there is no isopycnal T gradient, so maybe overinterpreting?
TODO:
[x] compute isopycnal gradients
[ ] Make sure I understand the negative sign on Kuz, Kvz
[ ] compare χ from Cole to ECCO K adjusted map
[ ] I have trouble matching the published, remapped Kredi to Kux; issues might be the region I picked. Should get the data on the native grid.
kredi.sel(k=1000, method="nearest").sel(lat=slice(20, 40), lon=slice(290, 360)).hvplot()
(ecco.GM_Kux).cf.plot.line(hue="node", y="Z")
(ecco.GM_Kvy).cf.plot.line(hue="node", y="Z")
ecco.Ke.cf.plot(y="k", color="k")
[<matplotlib.lines.Line2D at 0x1783ec820>]
ecco.GM_CHI.cf.plot(hue="node", y="vertical")
[<matplotlib.lines.Line2D at 0x177eb07c0>,
<matplotlib.lines.Line2D at 0x177f43a60>]
3-term balance?#
pop_var = xr.open_zarr("../datasets/pop-variance-budget.zarr/")
pop_var
<xarray.Dataset>
Dimensions: (depth: 62, lat: 2400, lon: 3600)
Coordinates:
TLAT (lat, lon) float64 dask.array<chunksize=(400, 800), meta=np.ndarray>
TLONG (lat, lon) float64 dask.array<chunksize=(400, 800), meta=np.ndarray>
* depth (depth) float64 0.05 0.15 0.25 0.35 ... 51.25 53.75 56.25 58.75
Dimensions without coordinates: lat, lon
Data variables:
BC (depth, lat, lon) float64 dask.array<chunksize=(62, 400, 800), meta=np.ndarray>
HDIFF (depth, lat, lon) float64 dask.array<chunksize=(62, 400, 800), meta=np.ndarray>
PKC (depth, lat, lon) float64 dask.array<chunksize=(62, 400, 800), meta=np.ndarray>
VMIX (depth, lat, lon) float64 dask.array<chunksize=(62, 400, 800), meta=np.ndarray>SSH urms / Salinity on isopycnals#
argoclim = xr.open_zarr(
"../datasets/argo-cole-clim-sigma2-2023-04-21.zarr/"
).cf.guess_coord_axis()
argoclim
<xarray.Dataset>
Dimensions: (lat: 145, lon: 360, σ: 60)
Coordinates:
* lat (lat) float64 -64.5 -63.5 -62.5 -61.5 -60.5 ... 76.5 77.5 78.5 79.5
* lon (lon) float64 20.5 21.5 22.5 23.5 24.5 ... 376.5 377.5 378.5 379.5
* σ (σ) float32 34.15 34.15 34.17 34.18 34.22 ... 37.2 37.2 37.2 37.2
Data variables:
CT (lat, lon, σ) float64 dask.array<chunksize=(50, 50, 60), meta=np.ndarray>
SA (lat, lon, σ) float64 dask.array<chunksize=(50, 50, 60), meta=np.ndarray>
z_σ (lat, lon, σ) float32 dask.array<chunksize=(50, 50, 60), meta=np.ndarray>hires = xr.open_zarr(
"../datasets/pop-hires-annual-climatology-grid-vars.zarr",
chunks={"nlat_t": 600, "nlon_t": 600},
)
hires["TLONG"] = xr.where(hires.TLONG > 180, hires.TLONG - 360, hires.TLONG)
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/dataset.py:255: UserWarning: The specified Dask chunks separate the stored chunks along dimension "nlon_t" starting at index 600. This could degrade performance. Instead, consider rechunking after loading.
warnings.warn(
import pop_tools
lores = xr.load_dataset("../datasets/pop-1deg-sigma2-42-61-clim.nc")
grid, logrid = pop_tools.to_xgcm_grid_dataset(pop_tools.get_grid("POP_gx1v7"))
lores["TLAT"] = grid.interp(logrid.ULAT, ["X", "Y"])
lores["TLONG"] = grid.interp(logrid.ULONG, ["X", "Y"])
lores["TLONG"] = xr.where(lores.TLONG > 180, lores.TLONG - 360, lores.TLONG)
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import colorcet
kwargs = dict(
y="latitude",
x="longitude",
vmin=35,
vmax=36,
add_colorbar=False,
)
f, ax = plt.subplots(
1,
4,
constrained_layout=True,
subplot_kw=dict(projection=ccrs.Orthographic(-20, 20)),
)
def plot_salt(da, ax):
h = da.isel(σ=36).cf.plot(ax=ax, cmap="cet_bmy_r", **kwargs)
da.isel(σ=36).cf.plot.contour(
ax=ax, colors="darkgray", linewidths=1, levels=np.linspace(35, 36, 11), **kwargs
)
return h
plot_salt(argoclim.SA.sel(lat=slice(10, 50), lon=slice(360 - 40, 360)), ax[0])
plot_salt(hires.SALT.isel(nlat_t=slice(1250, 1650), nlon_t=slice(700, 1140)), ax[1])
plot_salt(lores.SALT.cf.isel(X=slice(0, 35), Y=slice(225, 320)), ax[2])
h = plot_salt(ecco_tile2.SALT, ax=ax[3])
[ed.plot.make_nice_natre_map(ax_) for ax_ in ax]
ax[0].set_title("(a) Argo")
ax[0].gridlines(draw_labels=["left", "bottom"])
ax[1].set_title("(b) POP 1/10°")
ax[1].gridlines(draw_labels=["bottom"])
ax[2].set_title("(c) POP 1°")
ax[2].gridlines(draw_labels=["bottom"])
ax[3].set_title("(d) ECCO")
ax[3].gridlines(draw_labels=["right", "bottom"])
f.colorbar(
h,
ax=ax,
orientation="horizontal",
extend="both",
label="Salinity on $σ_{2}$=36.33 kg/m3",
aspect=40,
shrink=0.6,
)
f.set_size_inches((10, 4))
f.savefig("../images/isopycnal-salt.png", dpi=200, bbox_inches="tight")
2023-04-26 22:15:52,337 - tornado.application - ERROR - Exception in callback <bound method SystemMonitor.update of <SystemMonitor: cpu: 1 memory: 15 MB fds: 33>>
Traceback (most recent call last):
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/tornado/ioloop.py", line 921, in _run
val = self.callback()
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/system_monitor.py", line 161, in update
net_ioc = psutil.net_io_counters()
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/psutil/__init__.py", line 2114, in net_io_counters
rawdict = _psplatform.net_io_counters()
OSError: [Errno 12] Cannot allocate memory
2023-04-27 08:50:02,100 - distributed.scheduler - WARNING - Worker failed to heartbeat within 300 seconds. Closing: <WorkerState 'tcp://127.0.0.1:59063', name: 2, status: running, memory: 0, processing: 0>
2023-04-27 08:50:02,330 - distributed.scheduler - WARNING - Worker failed to heartbeat within 300 seconds. Closing: <WorkerState 'tcp://127.0.0.1:59065', name: 0, status: running, memory: 0, processing: 0>
2023-04-27 08:50:02,417 - distributed.scheduler - WARNING - Received heartbeat from unregistered worker 'tcp://127.0.0.1:59065'.
2023-04-27 08:50:02,429 - distributed.scheduler - WARNING - Received heartbeat from unregistered worker 'tcp://127.0.0.1:59063'.
2023-04-27 08:50:02,436 - distributed.worker - ERROR - Scheduler was unaware of this worker 'tcp://127.0.0.1:59063'. Shutting down.
2023-04-27 08:50:02,434 - distributed.worker - ERROR - Scheduler was unaware of this worker 'tcp://127.0.0.1:59065'. Shutting down.
Traceback (most recent call last):
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/utils.py", line 752, in wrapper
return await func(*args, **kwargs)
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/worker.py", line 1525, in close
await self.finished()
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/core.py", line 502, in finished
await self._event_finished.wait()
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/locks.py", line 226, in wait
await fut
asyncio.exceptions.CancelledError
Traceback (most recent call last):
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/utils.py", line 752, in wrapper
return await func(*args, **kwargs)
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/worker.py", line 1525, in close
await self.finished()
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/core.py", line 502, in finished
await self._event_finished.wait()
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/locks.py", line 226, in wait
await fut
asyncio.exceptions.CancelledError
2023-04-27 08:50:06,666 - distributed.nanny - ERROR - Worker process died unexpectedly
Exception in thread Nanny stop queue watch:
Traceback (most recent call last):
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/threading.py", line 973, in _bootstrap_inner
2023-04-27 08:50:06,679 - distributed.nanny - ERROR - Worker process died unexpectedly
2023-04-27 08:50:08,541 - distributed.nanny - WARNING - Restarting worker
2023-04-27 08:50:08,693 - distributed.nanny - WARNING - Restarting worker
2023-04-27 08:50:29,898 - distributed.nanny - WARNING - Worker process still alive after 3.1999897766113286 seconds, killing
2023-04-27 08:50:29,963 - distributed.nanny - WARNING - Worker process still alive after 3.199994201660157 seconds, killing
2023-04-27 08:50:29,998 - distributed.nanny - ERROR - Error in Nanny killing Worker subprocess
Traceback (most recent call last):
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/tasks.py", line 452, in wait_for
fut.result()
asyncio.exceptions.CancelledError
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 607, in close
await self.kill(timeout=timeout, reason=reason)
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 390, in kill
await self.process.kill(reason=reason, timeout=0.8 * (deadline - time()))
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 844, in kill
await process.join(max(0, deadline - time()))
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/process.py", line 330, in join
await wait_for(asyncio.shield(self._exit_future), timeout)
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/utils.py", line 1849, in wait_for
return await asyncio.wait_for(fut, timeout)
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/tasks.py", line 454, in wait_for
raise exceptions.TimeoutError() from exc
asyncio.exceptions.TimeoutError
2023-04-27 08:50:30,100 - distributed.nanny - ERROR - Error in Nanny killing Worker subprocess
Traceback (most recent call last):
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/tasks.py", line 452, in wait_for
fut.result()
asyncio.exceptions.CancelledError
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 607, in close
await self.kill(timeout=timeout, reason=reason)
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 390, in kill
await self.process.kill(reason=reason, timeout=0.8 * (deadline - time()))
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 844, in kill
await process.join(max(0, deadline - time()))
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/process.py", line 330, in join
await wait_for(asyncio.shield(self._exit_future), timeout)
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/utils.py", line 1849, in wait_for
return await asyncio.wait_for(fut, timeout)
File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/tasks.py", line 454, in wait_for
raise exceptions.TimeoutError() from exc
asyncio.exceptions.TimeoutError
aviso = xr.open_zarr("~/datasets/aviso_monthly.zarr", consolidated=True)
aviso["longitude"] = aviso.longitude - 360
aviso = aviso.sel(latitude=slice(10, 50), longitude=slice(-40, 5))
aviso["urms"] = np.sqrt(aviso.ugosa**2 + aviso.vgosa**2)
aviso["urms_mean"] = aviso.urms.mean("time").load()
aviso.urms_mean
<xarray.Dataset>
Dimensions: (time: 312, latitude: 160, longitude: 160, nv: 2)
Coordinates:
* latitude (latitude) float32 10.12 10.38 10.62 10.88 ... 49.38 49.62 49.88
* longitude (longitude) float32 -39.88 -39.62 -39.38 ... -0.625 -0.375 -0.125
* nv (nv) int32 0 1
* time (time) datetime64[ns] 1993-01-01 1993-02-01 ... 2018-12-16
Data variables:
adt (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
crs int32 ...
lat_bnds (latitude, nv) float32 dask.array<chunksize=(160, 2), meta=np.ndarray>
lon_bnds (longitude, nv) float32 dask.array<chunksize=(160, 2), meta=np.ndarray>
sla (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
ugos (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
ugosa (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
vgos (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
vgosa (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
urms (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
urms_mean (latitude, longitude) float64 0.1002 0.1002 ... 0.1431 0.1444
Attributes: (12/45)
Conventions: CF-1.6
Metadata_Conventions: Unidata Dataset Discovery v1.0
NCO: netCDF Operators version 4.7.9 (Homepage...
cdm_data_type: Grid
comment: Sea Surface Height measured by Altimetry...
contact: servicedesk.cmems@mercator-ocean.eu
... ...
summary: SSALTO/DUACS Delayed-Time Level-4 sea su...
time_coverage_duration: P1D
time_coverage_end: 1993-01-01T12:00:00Z
time_coverage_resolution: P1D
time_coverage_start: 1992-12-31T12:00:00Z
title: DT merged all satellites Global Ocean Gr...import cartopy.crs as ccrs
f, ax = plt.subplots(
1,
2,
constrained_layout=True,
subplot_kw=dict(projection=ccrs.Orthographic(-20, 20)),
)
aviso.urms_mean.plot(
transform=ccrs.PlateCarree(),
robust=True,
cmap="Reds",
vmin=0,
ax=ax[0],
add_colorbar=False,
vmax=0.15,
)
ed.plot.make_nice_natre_map(ax[0])
ax[0].set_title("(a) AVISO")
ax[0].gridlines(draw_labels=["left", "bottom"])
<cartopy.mpl.gridliner.Gridliner at 0x16df1a790>
Vertical grid spacings#
ECCO#
ecco.drC.set_xindex("Zp1").sel(Zp1=slice(800, 1500))
<xarray.DataArray 'drC' (k_p1: 7)>
array([ 95.27 , 97.415, 98.75 , 99.63 , 100.67 , 102.945, 107.945],
dtype=float32)
Coordinates:
* k_p1 (k_p1) int64 27 28 29 30 31 32 33
face int64 2
* Zp1 (k_p1) float32 861.5 958.0 1.056e+03 ... 1.357e+03 1.461e+03
drC (k_p1) float32 95.27 97.42 98.75 99.63 100.7 102.9 107.9
PHrefF (k_p1) float32 8.451e+03 9.398e+03 ... 1.331e+04 1.434e+04
Attributes:
standard_name: cell_z_size_at_w_location
long_name: cell z size
units: mPOP 1/10°#
pop_natre["hires/clim"].sel(depth=slice(800, 1500))["depth"].diff("depth")
<xarray.DataArray 'depth' (depth: 4)>
array([105.88061523, 121.49841309, 138.3626709 , 155.93029785])
Coordinates:
* depth (depth) float64 984.7 1.106e+03 1.245e+03 1.4e+03
time object 0052-01-01 05:17:48.750000
σ (depth) float64 36.57 36.68 36.77 36.86
z_σ (depth) float64 984.7 1.106e+03 1.245e+03 1.4e+03
Attributes:
positive: down
standard_name: depth
units: meterPOP 1°#
pop_natre["1deg"]
<xarray.DatasetView>
Dimensions: ()
Data variables:
*empty*Error bounds#
Microstructure#
Ferrari & Polzin (2005):

Studying Cole estimate#
Previously I was seeing some big disagreement. This was due to gradients estimated in depth-space instead of in isopycnal space.
Compare Cole & Groeskamp diffusivities#
Reproducing Groeskamp Figure 1. This seems to work well; but note the different domain. I saw this in his MATLAB script.
%% MIXING IN NATRE REGION
% THE NATRE REGIONS LIES BETWEEN
% 319.5 - 344.5 degrees lon
% 19.5 - 29.5 degrees lat
NATRE_x = 320:345;
NATRE_y = 110:120;
cole = ed.read_cole(resolution="3deg").sel(
lat=slice(19.5, 29.5), lon=slice(319.5, 344.5)
)
groeskamp = (
xr.open_dataset("/home/deepak/datasets/groeskamp2020/groeskamp2020.nc")
.sel(lat=slice(19.5, 29.5), lon=slice(319.5, 344.5))
.rename({"depth": "pres"})
)
def plot_mean_std(da, label):
mean = da.mean(["lat", "lon"])
std = da.std(["lat", "lon"])
(hdl,) = mean.cf.plot(ylim=(2000, 0), label=label)
dcpy.plots.fill_between(
xr.concat((mean - std, mean + std), dim="err"),
axis="x",
x="err",
y="pres",
alpha=0.3,
)
plot_mean_std(groeskamp.Ke_0, label="Groeskamp unsuppressed")
plot_mean_std(groeskamp.Ke, label="Groeskamp suppressed")
plot_mean_std(cole.diffusivity, label="Cole")
plt.gcf().set_size_inches((3.5, 5))
plt.gca().set_xlim((-500, 3000))
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5))
/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice.
var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
<matplotlib.legend.Legend at 0x7f6a3f8c57f0>
Depth-vs-density space estimation of \(K_e |∇T|²\)#
Turns out this does not matter; I was worried that diffusivity in depth space and |∇T| estimated in density space might combine badly. But moving diffusivity to sigma space or |∇T| to depth space results in same mean variance tendency
cole = ed.read_cole()
import xgcm
cole["diffusivity_sig"] = xgcm.transform.linear_interpolation(
cole.diffusivity, cole.density_mean_depth, cole.sigma, "pres", "pres", "sigma"
)
# read, take mean in region; swap to pressure coordinate
colegrad = xr.open_dataset("../datasets/cole-clim-gradient-natre.nc").sel(
lat=slice(24, 28), lon=slice(360 - 30, 360 - 25)
)
colegrad["sigma"] = cole.sigma
cole_natre = cole[
["diffusivity", "diffusivity_sig", "density_mean_depth", "depth_mean_sig"]
].interp(lat=colegrad.lat.data, lon=colegrad.lon.data)
dTiso_pres = xgcm.transform.linear_interpolation(
colegrad.dTiso,
colegrad.pres,
cole.pres,
"sigma",
"sigma",
"pres",
)
cole_var_sig = (
(cole_natre.diffusivity_sig * colegrad.dTiso**2)
.mean(["lat", "lon"])
.assign_coords(pres=cole_natre.depth_mean_sig.mean(["lat", "lon"]))
)
cole_var_pres = (cole_natre.diffusivity * dTiso_pres**2).mean(["lat", "lon"])
cole_var_sig.plot(y="pres")
cole_var_pres.plot(y="pres", yincrease=False)
plt.legend(["sigma space", "depth space"])
<matplotlib.legend.Legend at 0x7fa53eb855b0>
Depth space calculation#
The following is when I was using \(|∇T|\) estimated in depth space and multiplying by a \(K_e\) provided in depth space
Initial surprising observation: We expect to see peak \(K_e\) between 600 and 1600m but we mostly don’t! Both \(⟨S'S'⟩\) and \(|∇S|\) are large in the depth we expect a peak. There seems to be approximate agreement in the 1000-1200m range at (-34°E, -32°E) but missing data at deeper depths
Here are estimates of variance produced by eddy stirring at points in the NATRE box from Cole et al (2015). There is some spread but the peak between 400 and 1000 seems robust. Thick black line is the mean over the box.
cole_var_all.stack(latlon=["lat", "lon"]).drop("latlon").cf.plot.line(
y="pres", add_legend=False, size=4, aspect=0.8, xscale="log"
)
cole_var.cf.plot(color="k", lw=2)
plt.xlabel("$K_e^{cole} |∇T_{iso}^{argo}|²$ [°C²/s]")
# (cole.diffusivity.interp_like(dTiso) * colegrad.dSiso**2).mean(["lat", "lon"]).cf.plot()
Text(0.5, 0, '$K_e^{cole} |∇T_{iso}^{argo}|²$ [°C²/s]')
These are the various terms in the Cole analysis. Note minimum in ∇S between 200 and 800, and associated high mixing lengths ≳ 300km.
ed.plot.plot_cole_profile(lat=26, lon=360 + np.array([-35, -34, -33, -32, -30, -27]))
T-S diagrams: Argo vs. NATRE#
Lets compare T-S to see if they are similar. I am exploring ways to
interpret the “high” variance production values between 400m and 800m for the Cole estimate. Maybe the NATRE T-S is not representative of the argo T-S?
It is representative!
explain why the Argo ⟨χ⟩ seems to miss any eddy contribution. Is the scatter less at those density levels?
There is some undersampling (of course) because the sampling dz is low below 1000m usually.
It seems to capture some of the spread; but the NATRE dataset seems to fill up the space more.
The density contours are ρ corresponding to depths np.arange(200, 2000, 200) in the Cole dataset. Clearly the top few contours see not much T-S stirring and the T-S relationship is “tight”. Ferrari & Polzin interpreted this as a balance between microscale stirring and dissipation.
So the high K_e estimate between 400m and 800m in Cole et al (2015) is not supported by the TS spread.
Though interestingly the Groeskamp estimate has the same issue then.
Is the finestructure estimate undersampling eddy stirring periods? Compare T-S spread for
all argo profiles
argo profiles where I can make a finestructure estimate.
NATRE profiles
Show code cell source
argo_natre = xr.load_dataset("../datasets/argo/natre.nc")
argo_natre["THETA"] = dcpy.eos.ptmp(
argo_natre.PSAL, argo_natre.TEMP, argo_natre.PRES, natre.reference_pressure.item()
)
argo_fine_full_profiles = xr.load_dataset("argo_profiles_used_for_finestructure.nc")
# choose density levels that highlight the depth range we are interested in
rho_levels = (
cole.density_mean_depth.interp_like(colegrad)
.mean(["lat", "lon"])
.sel(pres=np.arange(200, 2000, 200), method="nearest")
)
All argo vs NATRE#
Show code cell source
kwargs = dict(
hexbin=False,
size=0.5,
plot_kwargs={"alpha": 0.1},
Sbins=np.arange(35, 36, 0.05),
Tbins=np.arange(2, 15, 0.5),
Pref=0,
)
_, ax = dcpy.oceans.TSplot(
argo_natre.PSAL,
argo_natre.THETA,
**kwargs,
rho_levels=rho_levels,
)
# mask = argo_fine_full_profiles.fine_deepest_bin > 1000
# _, ax = dcpy.oceans.TSplot(
# argo_fine_full_profiles.where(mask, drop=True).PSAL,
# argo_fine_full_profiles.where(mask, drop=True).THETA,
# **kwargs,
# rho_levels=rho_levels,
# )
dcpy.oceans.TSplot(
natre.salt,
natre.theta,
ax=ax,
**kwargs,
)
ax["s"].legend(["argo full", "natre"])
plt.gcf().savefig("../images/ts-natre.png")
Finestructure Argo vs NATRE#
Show code cell source
kwargs = dict(
hexbin=False,
size=0.5,
plot_kwargs={"alpha": 0.1},
Sbins=np.arange(35, 36, 0.05),
Tbins=np.arange(2, 10, 0.5),
Pref=0,
)
The distributions look similar and the argo data capture more extremes; but the NATRE data fill a lot of the space (though smaller S.D.). Maybe the average isn’t that great and it’s totally undersampling
Show code cell source
mask = argo_fine_full_profiles.fine_deepest_bin > 1000
_, ax = dcpy.oceans.TSplot(
natre.salt,
natre.theta,
**kwargs,
rho_levels=rho_levels,
)
dcpy.oceans.TSplot(
argo_fine_full_profiles.where(mask, drop=True).PSAL,
argo_fine_full_profiles.where(mask, drop=True).THETA,
ax=ax,
**kwargs,
)
ax["s"].legend(["natre", "argo fine"])
plt.gcf().savefig("../images/ts-natre.png")
It’s definitely a big undersampling when compared to all Argo profiles in the area. How good is good enough?
Show code cell source
_, ax = dcpy.oceans.TSplot(
argo_natre.PSAL,
argo_natre.THETA,
**kwargs,
rho_levels=rho_levels,
)
mask = argo_fine_full_profiles.fine_deepest_bin > 1000
masked = argo_fine_full_profiles.where(mask, drop=True)
_, ax = dcpy.oceans.TSplot(
masked.PSAL,
masked.THETA,
ax=ax,
rho_levels=None,
**kwargs,
)
masked.groupby
ax["s"].legend(["argo full", "argo fine"])
plt.gcf().savefig("../images/ts-natre-argo.png")
Conclusions#
Can we define a metric using (T-S spread in all argo profiles) vs (TS spread in profiles where we can make a finestructure estimate) to decide where we might be able to make an estimate using Argo?
Can we refine this?
If we expect regions above 1000m then it might not be so bad.
Could make a bigger box? (skeptical)
This really points out the limitation of the finestructure estimate, at least deeper down?
How does the Whalen scatter plot work out?
Is there spatial uncertainty? Cole estimate: Are we looking in the wrong place? No#
Here’s the NATRE region marked on a map of \(K_e\). There is just low diffusivity everywhere in the box; but may be the mean signal is just northwest of the box
fg = (
cole.diffusivity.sel(pres=[800, 1200], method="nearest")
.sel(lon=slice(300, 350), lat=slice(20, 40))
.plot(robust=True, row="pres", cbar_kwargs={"orientation": "horizontal"})
)
fg.map(
lambda: plt.plot(
360 + np.array([-31, -26, -26, -31, -31]), [24, 24, 28, 28, 24], lw=2
)
)
<xarray.plot.facetgrid.FacetGrid at 0x7fbb1793f2b0>